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ABSTRACT 

We present new measurements of the temperature of the intergalactic medium 
(IGM) derived from the Lya forest over 2.0 < z < 4.8. The small-scale structure in 
the forest of 61 high-resolution QSO spectra is quantified using a new statistic, the 
curvature, and the conversion to temperature calibrated using a suite of hydrodynamic 
simulations. At each redshift we focus on obtaining the temperature at an optimal 
overdensity probed by the Lya forest, T(A), where the temperature is nearly a one- 
to-one function of the curvature regardless of the slope of the temperature-density 
relation. The median 2er statistical uncertainty in these measurements is 8 per cent, 
though there may be comparable systematic errors due to the unknown amount of 
Jeans smoothing in the IGM. We use our T(A) results to infer the temperature at 
the mean density, To. Even for a maximally steep temperature-density relation, To 
must increase from ~8 000K at z ~ 4.4 to >12 000K at z ~ 2.8. This increase is 
not consistent with the monotonic decline in To expected in the absence of He n 
reionization. We therefore interpret the observed rise in temperature as evidence of 
He II reionization beginning at z > 4.4. The evolution of To is consistent with an end 
to He II reionization at z ~ 3, as suggested by opacity measurements of the He n 
Lya forest, although the redshift at which To peaks will depend somewhat on the 
evolution of the temperature-density relation. These new temperature measurements 
suggest that the heat input due to the reionization of He n dominates the thermal 
balance of the IGM over an extended period with Az > 1. 

Key words: intergalactic medium - quasars: absorption lines - cosmology: observa- 
tions - cosmology: early Universe 



1 INTRODUCTION 

Cosmic reionization represents one of the most dramatic 
feedback processes occurring between luminous sources and 
the intergalactic medium (IGM). At some point before red- 
shift six, nearly every hydrogen atom in the Universe was 
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ionized, as evidenced by the tr ansmission of H I Lya photons 
in Q S O spe ctra up to z ~ 6 (|Becker et al.l [200ll : iFan et all 



1200 j . bOOrj ). Similarly, helium in the IGM is known to be 
fully ionized by z ~ 2.7 based on observed transmission in 
the He II L ya forest at z < 3 (for a summary of obser- 



Dixon fc Furlanettd[2009l 'l. Understanding how 



vations, see 

and when these reionization events occurred is fundamen- 
tal to understanding the nature of the earliest galaxies and 
AGN, as well in assessing the impact of reionization on the 
evolution of luminous sources up to the present day. 

One of the main impacts of hydrogen and helium 
reionization is on the thermal state of the IGM. This is 
particularly true at low densities, where photoionization 
heating is believed to be the dominant heating process. 
When the bulk of the hydrogen and helium in the IGM 
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becomes ionized, the heat input is expected to signifi- 
cantly increase the gas temperature. The thermal history 
of the IGM will reflect not only the timing and duration 
of reionization, but also the nature of the ionizing sources 
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2009). At the mean den- 



sity, the characteristic signature of reionization is expected 
to be a marked heating followed by a relatively rapid cooling, 
the latter du e to adiabatic expansion and (a t z > 6) Comp- 
ton c ooling (|Miralda-Escude fc Reesl 1 1994 iHui fc Gnedinl 
Il997ft . If helium reionizes later than hydrogen, as expected 
if helium is fully reionized by bright QSOs, then the temper- 
ature evolution will contain two peaks. Following reioniza- 
tion, the thermal evolution will asymptotically approach a 
gradual cooling state, in which adiabatic expansion is nearly 
balanc ed by photoionization heating from H I, H e I, and/or 
He II (|Hui fc Gnedinl Il997l ; IHui fc Haimanll2003l ). In short, 
the thermal history of the IGM reflects both the output of 
luminous sources and the growth of large-scale structure. 

At z < 6, the most readily observable impact of the 
IGM temperature is on the small-scale structure of the 
H I Lya forest. The widths of absorption features are 
set by thermal broadening, Hubble broadening due to the 
change in redshift across a filament or void, peculiar veloc- 
ities, and the column densities of the absorbers. By using 
numerical simulations to capture the large-scale structure 
and bulk motions of the IGM, therefore, the temperature 
can be extracted by analyzing the amount of small-scale 
structure in the forest. Previous efforts to to do this have 
taken two main approaches. The first has been to mea- 
sure the distributio n of Doppler parameters of individua l 
absorption features jSchave et all l2000l ; iRicotti et al.l |2000| ; 
iMcDonald et alj200lft . The widths of the narrowest features 
are expected to be dominated by thermal broadening, and so 
can be used to infer a temperature. Alternatively, the small- 
scale structure has been quantified acro ss the entire for- 
est u sing either the flux power sp ectrum dZaldarriaga et al.l 
200 if) or wavelet decomposition (ITheuns fc Zaroubil l2000t 



Zaldarriagall2002l ; iTheuns et alj|2002al lbl). 

Obtaining reliable measurements of the IGM tem- 
perature, however, has remained a considerable challenge. 
ISchave et all (|2000h found tentative evidence for an increase 
in the IGM temperature consistent with He H r eionization 
near z ~ 3. In constrast, IMcDonald et al] (|200lh found the 
IGM temperature to be roughly constant over z ~ 2 — 4 us- 
ing a very si milar data set but a somewhat different line fit- 
ting method. iLidz et al] l|2009h recently conducted the most 
careful analysis to date using wavelets. Their study included 
a large data set and paid close attention to systematic ef- 
fects, yet the uncertainties in the temperature measurements 
remained substantial. 

This work aims to produce robust IGM temperature 
measurements over 2 < z < 5 using a method that is new in 
two ways. First, we introduce a new statistic, the curvature, 
that is well suited to quantifying the small-scale structure 
of the Lya forest. The curvature is essentially the second 
derivative of the flux with respect to wavelength or relative 
velocity, and has a number of advantages for this type of 
work. It is straightforward to compute and can be evaluated 



on a pixel-by-pixel basis. It also does not require the forest to 
be decomposed into individual lines, and so can be applied 
at z > 4 where absorption features start to become strongly 
blended. Second, we focus our efforts on accurately measur- 
ing the temperature at a well-chosen overdensity, rather than 
attempting to constrain the entire temperature-density rela- 
tion. The temperature-density relation is typically described 
by a power law, T(A) = To A 7-1 , where To is the tempera- 
ture at the mean density and A = p/(p). We show that, at 
each redshift, an optimal overdensity, A, can be identified 
where T(A) is a one-to-one function of the mean absolute 
curvature, regardless of 7. This allows us to normalize the 
temperature-density relation with high precision, providing 
a first step towards a full characterization of the thermal 
state of the IGM. 

The remainder of the paper is organized as follows: 
We introduce the data in Section [2] and the simulations to 
which the data are compared in Section [3] We introduce the 
curvature statistic and describe our method for translating 
the curvature into a temperature measurement in Section [4] 
Here we also discuss our strategies for dealing with sources 
of potential systematic error, as well as examine the poten- 
tial impact of the thermal history of the IGM on measur- 
ing instantaneous temperatures. Our results are presented 
in Section [5] These include the temperature at the optimal 
overdensity, as well as the temperature at the mean density 
for different values of 7. The implications for He 11 reion- 
ization are discussed in Section [6] Finally, we summarize 
our conclusions in Section [7] The numerical convergence of 
our simulations is examined in the Appendix. Comoving dis- 
tances are used throughout the paper. 



2 THE DATA 

Our sample includes 61 QSOs spanning emission redshifts 
2.1 < Zem < 5.8. The objects are listed in Table [1] High- 
resolution spec tra were obtained primarily with the HIRES 
spectrograph (|Vogt et al.l Il994h on Keck I. Observations 
were made with the original CCD detector between 1993 
and 2004, and with the upgraded detector between 2005 
and 2008. Additional spectra of QSOs at z > 4.5 were taken 
with the MIKE spectrograph on Magellan (|Bernstein et al] 
12002ft between 2006 and 2008. For QSOs at 2 < 4.7 taken 
with HIRES, the raw data were processed and the extracted 
spectra combined using the MAKEE reduction package. At 
z > 4.7, the HIRES data wer e reduced using a cu stom set of 
IDL routines, as described in iBecker et al] (120061 ). The IDL 
routines were designed t o handle fain t targets, and include 
optimal sky subtraction (|Kelsonll2003ft as well as simultane- 
ous optimal extraction of the ID spectrum from all expo- 
sures. The MIKE data were reduced using a similar custom 
software package. 

The resolution of the HIRES data is 6.7 km s _1 
(FWHM). The majority of the HIRES data were binned 
using pixels with a constant velocity width of 2.1 km s _1 , 
although a subset of the spectra at z < 2.8 were binned using 
constant 0.03 A pixels. The resolution of the MIKE spectra 
is 13.6 km s~ . For these spectra, 5.0 km s _1 binned pixels 
were used. For both instruments the resolution is sufficient 
to resolve typical Lya forest lines, which have Doppler pa- 
rameters b > 10 km s~\ or FWHM > 24 km s -1 . When 
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Table 1. List of QSOs analyzed in this work. Columns give the 
QSO name, the emission redshift, the instrument with which the 
spectrum was taken, the total exposure time across all settings, 
and the mean noise in the normalized flux per 2.1km s _1 among 
those sections of the spectra that were included in the analysis. In 
order to be included, a 10 h^ 1 Mpc section was required to have 
a mean noise level less than 0.06. Objects for which 'HIRES' 
is marked with an (*) were observed at least partially with the 
upgraded detector. 

QSO z em Instrument t exp Noise 



Table 1. - Continued 



QSO 




Instrument 


£exp 


Noise 


SDSS 0040-0915 


4.98 


MIKE 


8.3 


0.054 


SDSS 0915+4924 


5.20 


HIRES* 


10.8 


0.058 


SDSS 1659+2709 


5.32 


HIRES* 


11.7 


0.041 


SDSS 0836+0054 


5.80 


HIRES* 


21.7 


0.058 



making comparisons with the simulations, we smooth and 
resample the simulations to match the resolution and bin- 
ning of each QSO spectrum individually. Provided that ad- 
equate sampling is used, we found the curvature statistic 
(Section [4} does not depend on the choice of binning, al- 
though it does depend on the resolution. As described below, 
the curvature is measured in sections spanning 10 ft -1 Mpc 
sections. Only sections with a mean rms uncertainty in 
the normalized flux of 0.06 per 2.1 km s" 1 birQ were in- 
cluded in the final analysis. Initial continua were fit using 
slowly-varying spline profiles. Real data are difficult to con- 
tinuum fit accurately, however, and we discuss our strategy 
for avoiding errors due to continuum fitting in Section ^. 2, 21 



3 HYDRO DYNAMICAL SIMULATIONS 

In order to obtain temperature constraints from the Lya for- 
est, we have used detailed simulated spectra for calibrating 
and interpreting our observational measurements. We per- 
formed a large set of high resolution, fully hydrodynamical 
simulations which span a wide variety of possible thermal 
histories. The simulations were run using the parallel Tree- 
SPH code GADGET-3, which is an updated versi on of the 
publicly available code GADGET-2 (|Springej|2005h . and are 
summarized in Table [2] 

The fiducial simulation volume is a 10 Mpc periodic 
box containing 2 x 512 gas and dark matter particles. This 
resolution is chosen specifically to re solve the Lya forest at 
high redshift (jBolton fc Beckerll2009h . The simulations were 
all started at z = 99, w ith initial conditions ge nerated using 
the transfer function o f lEisenstein fc Hull] 19991 ). The cosmo- 
logical parameters are fi m = 0.26, 0\ = 0.74, n h h 2 = 0.023, 
h = 0.72, CTg = 0.80, n s = 0.96, consistent with recent stud- 
ies o f the cosmic microw ave background (|Reichardt et al.l 
120091 : Ijarosik et ai1l2010f l. The IGM is assumed to be of 
primordia l composition with a h elium fraction by mass of 
Y = 0.24 ([Olive fc Skillmanll2004r ). The gravitational soften- 
ing length was set to l/30 th of the mean linear interparticle 
spacing and star formation was included using a simplified 
prescription which converts all gas particles with overden- 
sity A = p/(p) > 10 3 and temperature T < 10 5 K into 
collisionless stars. A total of 17 outputs at different red- 
shifts are obtained from each simulation, spanning the range 
1.831 < 2 < 6.010. 

The gas in the simulations is assumed to be optically 
thin and in ionization equilibrium with a spatially uniform 
ultraviolet background (UVB). The UVB corresponds to the 
galaxi es and quasars emission model of lHaardt fc MadatJ 
(|200ll ). Hydrogen is reionized at z = 9 and gas with A < 10 

1 This is a typical pixel size for the HIRES data. For MIKE data, 
the noise threshold per 5.0 km s -1 pixel was 0.06 X ^/2. 1/5.0. 
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subsequently follows a tight power-law temperature-density 
relation, T = Hp A 7-1 , where Tp is the temperature of the 
IGM at mean density (jHui fc Gnedinlll997l : IValageas et al.l 
120021 ). In o rder to explore a variety of thermal histories, we 
rescale the lHaardt fc Madaul (|200ll ) photo- heating rates by 
different constants in the models. In each simulation we as- 
sume €j = CA«e™ 01 , where e, HM01 are the lHaardt fc Madaul 
|200ll ) photo- heating rates for species i = [Hi, He 1, He 11] 
and £, £ are constants listed in Table Our fiducial mod- 
els (Af 5-GI 5) use £ = 0.0, which produces a temperature- 
density relation with 7 ~ 1.5 at 2 = 3 (Table [2j|. A selec- 
tion of other models include a density-dependent rescaling of 
the photo-heating rates to explore the effect of the relation- 
ship between temperature and gas density on our results. In 
two further models, T15slow and T15fast, we use a density- 
independent scaling which varies with redshift. These sim- 
ulations have integrated thermal histories which are signif- 
icantly different from that of the HM01 reference model. 
We use these models to test the impact of Jeans smooth- 
ing on our results in Section 14.51 Finally, we also perform 
four simulations, Rl, R2, R3 and R4, which use same ther- 
mal history as model C15 but have different box sizes and 
mass resolutions. These simulations are used to test the nu- 
merical convergence of our results, which is discussed in the 
appendix at the end of this paper. The values of To and 7 
in the hydrodynamical simulations at z = 3 are listed in 
columns 7 and 8 of Table [21 and the thermal histories as a 
function of redshift for all simulations except the T and R 
series are shown in Figure [1] 

We construct synthetic Lya spectra at each redshift 
using the output o f our hydrodynamical simulations (e.g., 
iTheuns et ai1ll99cf ). The spectra are then modified in post- 
processing to match the characteristics of the observational 
data as closely as possible. We now turn to discuss this as- 
pect of our simulations in more detail, and introduce our 
method for obtaining our IGM temperature constraints. 





Figure 1. Thermal properties of the the main simulations used 
in this work. Lines are color coded according to To . Line styles de- 
note different 7's - solid: 7 ~ 1.5, dotted: 7 ~ 1.3, dashed: 7 ~ 1.0, 
dot-dashed: 7 ~ 0.7. Top panel: Temperature at the mean density 
as a function of redshift. Letters correspond to the first letter of 
the simulation name in Table [2] Middle panel: Slope+1 of the 
temperature-density relation, where T(A) = ToA 7-1 . Bottom 
panel: The curvature statistic computed as described in Section^] 
Here, the curvature has been calculated directly from the noise- 
free spectra. 



4 METHOD 

4.1 The curvature statistic 

We introduce a new statistic, the curvature, to efficiently 
measure the amount of small-scale structure in the Lya for- 
est. The curvature, k, is conventionally defined as 

F" 

K = ry- . (1) 

[1 + (F") 2 1 1 

We take the derivative of the flux with respect to the veloc- 
ity separation between pixels, which is measured in km s -1 . 
In practice, the change in flux per km s is generally small, 
and so the denominator is always close to unity. The curva- 
ture is therefore essentially equal to the second derivative of 
the flux. 

The curvature has several advantages over other mea- 
surements of the small-scale structure. It is easy to compute, 
and can be evaluated on a pixel-by-pixel basis. It is also a 
continuous function, and does not require the forest to be 
deconstructed into individual absorption features. This al- 
lows the curvature to be easily calculated at high redshifts 
(z > 4), where high levels of absorption make it difficult to 
identify individual absorption lines. 



Figure [2] demonstrates the sensitivity of the curvature 
to the IGM temperature. We plot simulated spectra at two 
redshifts, z — 2.553 and z = 4.586, drawn from simulations 
B15 and D15, for which the temperature at the mean den- 
sity differs by roughly a factor of two but the slope of the 
temperature-density relation is the same. At both redshifts, 
the flux profile drawn from the colder simulation (B15) is 
significantly more highly peaked, with larger differences be- 
tween adjacent minima and maxima. This translates into 
significantly larger curvature values for the colder simula- 
tion, where the peaks in the curvature correspond to the in- 
flection points in the flux. The overall amount of small-scale 
structure in a section of forest is measured by computing 
the mean absolute curvature, (\k\) among pixels within a 
fixed range range in flux (see Section l4.2,2p . Simulations are 
then used to convert (|«|) into a temperature, as described 
in Section |4, 31 

4.2 Sources of Uncertainty 

4.2.1 Noise 

Several issues must be addressed before the curvature can be 
used as a reliable probe of the IGM temperature. The first 
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Table 2. Simulation runs used in this work. Columns give the simulation name, the box size, the total number of particles, the gas 
particle mass, the scaling factors for the UVB photo- heating rates (see text for details), the median volume- weighted gas temperature 
at the mean density at z = 3, the slope+1 of the temperature-density relation at z = 3, where T(A) = TqA 7 , and notes on the 



simulations. 


Model L [h 


_i Mpc] 


Particles 


J1J gas [ UJ 


c 

s 




T„ z = 3 fKl 




Notes 


A15 


10 


2 X 


512 3 


9.2 x 10 4 


0.30 


0.00 


5 100 


1.52 


Fiducial T-A relation 


B15 


10 


2 X 


512 3 


9.2 x 10 4 


0.80 


0.00 


9 600 


1.54 


Fiducial T-A relation 


C15 


10 


2 x 


512 3 


9.2 x 10 4 


1.45 


0.00 


14 000 


1.54 


Fiducial T-A relation 


CIO 


10 


2 x 


512 3 


9.2 x 10 4 


1.45 


-1.00 


13 700 


1.02 


"Isothermal" f 7 1 - /\ relation 


D15 


10 


2 x 


512 3 


9.2 x 10 4 


2.20 


0.00 


18 200 


1.55 


Fiducial T-A relation 


D13 


10 


2 x 


512 3 


9.2 x 10 4 


2.20 


-0.45 


18 100 


1.32 


Flattened T-A relation 


D10 


10 


2 x 


512 3 


9.2 x 10 4 


2.20 


-1.00 


18 000 


1.03 


"Isothermal" T-A relation 


D07 


10 


2 x 


512 3 


9.2 x 10 4 


2.20 


-1.60 


17900 


0.71 


Inverted T-A relation 


E15 


10 


2 x 


512 3 


9.2 x 10 4 


3.10 


0.00 


22 500 


1.55 


Fiducial T-A relation 


E10 


10 


2 x 


512 3 


9.2 x 10 4 


3.10 


-1.00 


22 200 


1.04 


"Isothermal" T-A relation 


F15 


10 


2 x 


512 3 


9.2 x 10 4 


4.20 


0.00 


27000 


1.55 


Fiducial T-A relation 


G15 


10 


2 x 


512 3 


9.2 x 10 4 


5.30 


0.00 


31000 


1.55 


Fiducial T-A relation 


T15slow 


10 


2 x 


512 3 


9.2 x 10 4 


Varied 


0.00 


18 200 


1.53 


Test run with gradual temperature increase 


T15fast 


10 


2 x 


512 3 


9.2 x 10 4 


Varied 


0.00 


18 600 


1.54 


Test run with rapid temperature increase 


Rl 


10 


2 x 


256 3 


7.4 x 10 s 


1.45 


0.00 


13 700 


1.02 


Mass resolution convergence 


R2 


10 


2 X 


128 3 


5.9 x 10 6 


1.45 


0.00 


13 700 


1.02 


Mass resolution/box size convergence 


R3 


20 


2 x 


256 3 


5.9 x 10 6 


1.45 


0.00 


14 000 


1.54 


Box size convergence 


R4 


40 


2 x 


512 3 


5.9 x 10 6 


1.45 


0.00 


14 000 


1.54 


Box size convergence 




-0.004 I . , . , . , . i . I . , . , J , . , . , , l 

200 400 600 800 200 400 600 800 1000 1200 

v (km/s) v (km/s) 



Figure 2. A demonstration of the curvature statistic. Top panels: 10 ft. -1 Mpc sections of simulated Lya forest drawn at z = 2.553 (left) 
and z = 4.586 (right). The sight lines plotted in blue are red drawn from simulations B15 and D15, respectively, which differ by roughly a 
factor of two in Tq but have similar 7's. Bottom panels: The curvature statistic computed from these spectra, as defined in equation ((TJ. 
At both redshifts, the overall amplitude of re is greater at the lower temperature. 



of these is noise in the flux spectra, which will dominate the 
curvature measurement if it is calculated directly from even 
very high-quality data. Rather than compute the curvature 
from the raw spectra, therefore, we first fit a smoothly vary- 
ing b-spline to the flux. The curvature is then computed 
from the fit. We use cubic polynomials for the piece-wise 
fits, such that the curvature is continuous across the break 
points. The b-spline if fit adaptively, with the break points 
initially separated by 50 km s _1 . The fit is then evaluated 
around each break point, and additional points are added 
where the fit is poor in a chi-squared sense. This process 
is repeated until the fit converges, or until the break point 



spacing reaches a minimum value (6.7 km s" 1 for the HIRES 
data, 13.6 km s _1 for the MIKE data). 

Figure [3] shows examples of the b-spline fits to the 
data. The fits are naturally more highly constrained for high 
signal-to-noise data. In order to compare with the simula- 
tions, therefore, we add the same level of the noise to the 
simulated spectra as we find in the data, and compute the 
b-splines using the same method. A significant benefit of our 
adaptive fitting scheme is that the curvature of the final fit 
is only moderately sensitive to the amount of noise in the 
spectrum. We estimate that our flux error arrays are gener- 
ally accurate to 10-20%. Tests using simulations indicated 
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Figure 3. Top panels: Examples of the Lya forest in QSO spectra used in this study, here at z = 2.470 taken with HIRES (left), and at 
z = 4.484 taken with MIKE (right). Each section spans 10 /i" 1 Mpc. The histograms show the normalized spectra. The continuous lines 
are b-spline fits. Bottom panels: The curvature computed from the b-spline fits. 



that uncertainties in the flux error at this level will trans- 
late into errors in the temperature measurements that are 
smaller than our statistical errors. We have also verified that 
the measured curvature is not highly sensitive to moderate 
correlations in the flux errors, and that a constant error is an 
adequate approximation in cases where the error scales with 
the flux, as expected for bright sources. We found that use- 
ful curvation informaiton could be extracted from data with 
errors in the normalized flux as large as ~0.1 per 2.1 km s _1 
bin. To be conservative, however, we have limited our anal- 
ysis to regions of the data where the mean error within a 
10 /i -1 Mpc section is less than 0.06 per 2.1km s _1 . This 
was done to maximize the contrast in curvature between 
simulation runs with different temperatures, which tends to 
diminish as sharp features become washed out by noise, and 
to make it easier to identify metal lines (see Section 14. 2. 4[1 , 



4-2.2 Continuum level 

The curvature as defined by Eq. ([TJ depends linearly on the 
amplitude of the flux, which in turn will depend on the con- 
tinuum level. Unfortunately, the high levels of absorption in 
the Lya forest make it difficult to reliably establish the con- 
tinuum in the real data, especially at z > 4, where even the 
voids may absorb > 10 per cent of the flux. This issue has re- 
cently been addressed in the context of measuring the mean 
transm itted flux in the Lya forest. IFaucher-Giguere et alj 
(|200&f ) used simulations to estimate their likely error in the 
continuum placement as a function of redshift, and then 
applied a correction to their measured mean flux. This ap- 
proach will work so long as the optical depth of the voids is 
consistent. We wish to allow for changes in the temperature- 
density relation, however, which can have a strong effect on 
the o pacity of the voids (|Becker et al.l 120071 : iBolton et alj 

l20oah . 



We therefore circumvent the continuum issue by "re- 
normalizing" both the real data and the simulations. For 
each 10 h~ Mpc section of forest (the size of our simula- 
tion box), we divide the flux by the maximum value of 
the b-spline fit in that interval. For example, the spectra 
in Figure [3] at z = 2.470 and 4.484 would be divided by 
1.00 and 0.87, respectively. The curvature is then computed 
from the re-normalized fit. This procedure forces us to sacri- 
fice a small amount of sensitivity. With the same mean flux 
level, hotter simulations result in spectra with lower over- 
all curvature but more absorption at the transmission peaks 
than colder simulations. Re-normalizing the flux therefore 
increases the curvature more in the hotter simulations, de- 
creasing the contrast in curvature between temperatures. 
The advantage, however, is that the continuum placement 
is handled objectively and consistently, removing a poten- 
tially significant sourc e of s ystematic error. For example, 
IFaucher-Giguere et al.l (|2008l ) estimate a continuum place- 
ment error of ~14 per cent at z — 4. If left uncorrected, this 
would translate into an error in the curvature of the same 
amount, which would in turn produce an error in tempera- 
ture of ~ 2 000 K. 

We measure the mean absolute curvature, (|«|), of pixels 
where the re-normalized b-spline fit falls in the range 0.1 < 
F R1 ° < 0.9, where the "R10" superscript denotes that the 
re-normalization has been performed in 10 ft -1 Mpc sections. 
The lower bound omits saturated pixels, which do not con- 
tain any information on the temperature. The upper bound 
is chosen largely to avoid uncertainties in the the mean abso- 
lute curvature due to uncertainties in the overall mean flux. 
This is particularly a concern at z ~ 2 — 3, where the flux 
probability distribution function starts to become sharply 
peaked near F = 1. Because the flux profile tends to be fairly 
flat near the continuum at these redshifts, the pixels there 
will tend to have lower than average absolute curvature val- 
ues. Increasing the mean flux, even within the best currently- 
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Figure 4. Top panel: The effective optical depth measured from 
our data after the transmitted flux has been re-normalized in 
10 h~ 1 Mpc sections (see Section 14. 2.21 1 and corrected for metal 
absorption. Vertical error bars are 2<r. The solid line is the fit 
used to scale the optical depths in our simulations. The dotted 
lines show the 2a scatter of the data points about this fit. Bottom 
panel: The conventional effective optical depth recovered from the 

simulations once r 1 ^ 10 has been scaled to match the observations. 

en 

The solid and dotted lines correspond to the same data as in the 
top panel, but here the flux has not been re-normalized. We use 
the range bounded by the dotted lines as our 2a uncertainty in r e ff 
when calculating uncertainties in the temperature. For clarity we 
show only r e ff from run C15, since the differences for other runs 
will be much smaller than the 2cr unc ertainties. For comparison, 
T e ff from lFa uchcr-Gig uere et al.l j200Sf) is plotted as a dashed line. 

measured uncertainties (e.g., iFaucher-Gieuere et al.ll2008h . 
significantly increases the number of pixels near the con- 
tinuum, and will therefore systematically decrease By 
not including pixels with F R1 ° > 0.9, we avoid most of this 
potential uncertainty. 

4.2.3 Mean flux 

Our choice to exclude pixels near the continuum minimizes 
the sensitivity of the curvature measurement to the overall 
mean flux level. Nevertheless, we have attempted to control 
for the mean flux in the simulations as carefully as pos- 
sible. We use the re-normalized flux spectra to measure a 
re-normalized effective optical depth, r^f 10 = — ln(i r ' R10 ). 
We then adjust the H I ionization rate in the simulations 
such that rjj 10 in the artificial spectra matches rjg 10 in the 
real data. This approach has the advantage over match- 
ing the true mean flux in that it does not depend on cor- 
rectly estimating the continuum in the real data. We also in- 
clude a correction for metal absorption. The flux decrement 
due to metals in the Lya forest was estimated by measur- 



ing the mean decrement between the N v A1239, 1243 and 
Si IV A1394, 1403 emission lines of the QSOs in our sam- 
ple. We found a nearly constant decrement of 3 per cent 
over 1.8 < 2 < 4.3. The decrement becomes more difficult 
to measure at higher redshifts as this rest wavelength in- 
terval begins to include regions of atmospheric absorption. 
However, the uncertainty in the temperature associated with 
uncertainty in the mean flux becomes small at z > 3 (see 
Section [5. 2) . We therefore adopt a constant 3 per cent cor- 
rection to the mean Lya flux due to metals for all redshifts. 

The redshift evolution of tJj 10 is shown in the top panel 
of Figure U The points show the mean values in the data, 
corrected for metal absorption, in redshift bins of Az — 0.4 
along with the 2a uncertainties. The solid line is the value 
of Ttfl° we set in our simulations, while the dotted lines 
show the 2a scatter of the data points about this fit. Once 
the simulations are scaled to produce the correct rjg 10 , the 
mean flux without re-normalizing can be measured from the 
artificial spectra. This is useful for comparing our mean flux 
constraints to results from the literature. We plot r e a in 
the bottom panel of Figure 3J where r e s = — \n{F). For 
clarity, we only show r e fi for the C15 simulation. The recov- 
ered mean flux will be slightly different for simulations with 
different temperatures; however, the spread in values is rel- 
atively small compared to t he uncertainty about the tren d 
shown. We also plot T e g from lFaucher-Giguere et al.l (|200ct ). 
which agrees well with our measurements, particularly at 
z < 3, where our method is most sensitive to the mean flux 
(see Section f5 . 2 f) . 

4.2.4 Metal lines 

Metal lines are a potentially serious source of systematic er- 
ror in any measure of the small-scale structure in the Lya 
forest. The individual components of metal absorption lines 
syste ms tend to be sig nificantly narrower than Lya lines 
(e.g., IWolfe et al.l l2005f ) . If regions of the forest contami- 
nated by metal lines are not properly masked, therefore, the 
temperature measurements will be significantly biased to- 
wards low values. 

We identify and mask metal lines using a two-step pro- 
cess. First, we identify as many absorption systems as pos- 
sible either from very strong Lya absorption (i.e., damped 
Lya systems and Lyman limit systems) or metal lines that 
fall redward of the Lya forest. The latter may include high- 
ionization species such as C IV, low-ionization species such 
as O 1, Si 11, C 11, and Fe 11, or any combination of lines that 
can be identified as belonging to the same system. We then 
mask all regions in the forest that could plausibly contain 
metal lines at the same redshifts as these systems. 

The next step is to search for any remaining unidentified 
metal lines in the forest. At high spectral resolution, metal 
absorption lines tend to stand out from the Lya forest by 
having one or more narrow (b < 10—15 km s _1 ) components, 
or very sharp edges when saturated. We therefore mask any 
lines that appear conspicuously narrow or sharp, as well as 
any doublets (e.g., C IV, Si IV, Mg 11) or other groups of lines 
that can be identified as belonging to a common system. 

Metal lines are a potentially greater contaminant at 
lower redshifts, where there are fewer Lya lines. However, 
this also makes the metal lines easier to identify and remove. 
While no metal rejection scheme can be perfect, we expect 
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that we have removed the large majority of metal lines that 
may significantly contaminate the curvature measurement. 
Since we are computing the mean absolute curvature across 
all pixels, the signal from any remaining unidentified metal 
lines should be weak in comparison to the signal from the 
Lya forest. For example, we have estimated the impact of 
unidentified C IV systems, which generally cannot be identi- 
fied using other lines redward of the forest. Random C IV sys- 
tems were added to the simulated spectra at a level consis- 
tent with the observed num ber density and column density 
distribution over z ~ 2 — 4 (|Songailal 200lh . Even if we as- 
sume that we have missed all systems with N < 10 14 cm" 2 , 
where the upper limit corresponds to a flux decrement near 
unity, and that these systems are all narrow (b = 10 km s _1 ), 
in run C15 at z ~ 4.6 would only increase by ~0.01 dex, 
which is equivalent to a decrease in To of ~400 K. At z ~ 3.0, 
the increase in (\k\) would be ~0.02 dex, corresponding to 
a decrease in To of ~700 K. We regard this as a conserva- 
tive upper limit, since at z < 3 the Lya forest is transpar- 
ent enough that strong doublets can often be identified and 
masked. We therefore expect the curvature measurement to 
be strongly dominated by the genuine Lya forest, and that 
the bias from unidentified metal lines will be smaller than 
our measurement errors (see Section [5]). 

4.3 Optimal overdensities 

A fundamental challenge to measuring IGM temperatures 
is that the temperature will vary with density. Following 
rcionization the gas in the IGM will be in photoionization 
equilibrium. At low densities (A < 10), the heating and 
cooling of the gas are dominated by the adiabatic contrac- 
tion and expansion of the filaments and voids due to struc- 
ture formation, and by the photoionization heating and re- 
combination cooling from H I, He I and/or He II. The bal- 
ance between the cooling and heating processes has been 
shown to produce a relatively tight, gradually steepening 
power law temperature-density relation of the form T(A) = 
Tp A 7-1 , where Tp is the temperature at the mean density 
(|Hui fc Gnedinlll997h . The slope is expected to asymptoti- 
cally reach 7 — 1 « 0.6, with the normalization set by the 
hardness of the ionizing background. 

Reionization will complicate this picture, however. The 
photoheating from reionization will be roughly independent 
of density. If reionization were homogeneous and instanta- 
neous, therefore, the temperature-density relation immedi- 
ately afterward would be nearly isothermal. In reality, reion- 
ization is likely to be an extended and inhomogeneous pro- 
cess. During He 11 reionization, photoheating within growing 
He ill bubbles around bright QSOs should partially flatten 
the volume-averaged temperature-density relation as well 
as introduce significant spatial fluctuations in temperature 
llGleser et all 120051; iFurlanetto fe Obll2008l ; iMcQuinn et all 
l2009l ; iBolton et alj 1200^ 7 At a minimum, therefore, mea- 
surements of the thermal state of the IGM must allow for 
changes in both the normalization and the slope of the 
temperature-density relation. 

Even constraining both Tp and 7 in a volume-averaged 
sense is challenging, however, since the dynamic range of 
the Lya forest is limited to densities that produce mea- 
surable but non-saturated absorption. Conventionally, the 
temperature-density relation is normalized at the mean den- 



sity (A = 1). From an interpretive standpoint, this is con- 
venient in that a parcel of gas near the mean density will 
tend to remain near the m ean density as the IGM evolves 
(e.g., iHui fc Gnedir]|l997h . Measuring To as a function of 
redshift can therefore provide a straightforward indication 
of how much heating or cooling is occurring for a relatively 
consistent population of baryons in the IGM. 

The difficulty with measuring To, however, is that the 
Lya forest does not always trace gas near the mean density. 
Ignoring the effects of peculiar velocities and thermal broad- 
ening, the Lya optical depth at a given overdensity will scale 
as r(A) oc ' H~ 1 (z) rim oc (1 + z )4.B r -i T -o.7 A a-o.r{ 7 -i) ) 
where H(z) is the Hubble parameter, nm is the neutral hy- 
drogen density, and F is the H I photoionization rate (e.g ., 
IWeinberg et al.1 1 19971 ; iMcDonald fc Miralda-Escudel 120011 ). 
The strong redshift dependence makes it extremely chal- 
lenging to probe a constant range of overdensities over all 
redshifts. In practice, the Lya forest will tend to be sensitive 
to gas that produces optical depths near unity, which occurs 
near the mean density at z ~ 4 — 5, and at mild overdensities 
at z ~ 2 — 3. At lower redshifts, therefore, the measurement 
of To will become increasingly degenerate with 7, since most 
of the signal will be coming from densities above the mean. 

The degeneracy of To and 7 in setting the amount of 
small-scale structure in the Lya forest is demonstrated in 
Figure [5] Here we show the same simulated sections of Lya 
forest as in Figure [5] but with simulation D15 replaced by 
D10, which has a nearly isothermal temperature-density re- 
lation (7 ~ 1). The gas at A > 1 is colder in this sim- 
ulation than it would be in the case of 7 ~ 1.5. Conse- 
quently, at z = 2.553 the Lya forest in B15 and D10 appears 
nearly identical, even though To differs by roughly a factor 
of two. The similarity is also seen in the curvature values. At 
z = 4.586 the fluctuations in the forest are still significantly 
smoothed out in run D10 as compared to run B15, since the 
forest at this redshift is more sensitive to gas near the mean 
density than at lower redshifts. 

In principle one can take advantage of the fact that 
lower-column density absorbers tend to arise from lower- 
density gas in order to constrain both Tp and 7 (e.g., 
ISchave et ail |2000| ; iTheuns fc Zaroubil I2OO0T ). In practice, 
however, this has proven to be extremely challenging, and 
the error bars on both quantities have remained large. At 
z > 4 the forest becomes increasingly difficult to separate 
into discrete absorbers, posing a further obstacle from spec- 
ifying the full temperature-density relation from the Lya 
forest alone. 

Rather than attempt to measure both the normalization 
and the slope of the temperature-density relation, therefore, 
we use the curvature to measure the temperature at an opti- 
mal overdensity probed by the Lya forest. The optimal over- 
density, A, which will change with redshift, is empirically 
determined to be the overdensity at which T(A) is a simple 
function of regardless of 7. Our approach is illustrated 
in Figure [6] In the left-hand panels we plot To as a function 
of Ioq(\k\) for the various simulation runs. The simulations 
with 7 ~ 1.5 form a well-defined relationship between the 
curvature and To. However, runs with different 7's can be 
well off of this relationship, particularly at lower redshifts. 
For a given log{\n\) at z ~ 2, the difference in To between 
7 ~ 1.5 and 7 ~ 1.0 can be > 10 000 K. This is clearly 
undesirable when 7 is poorly constrained. In the right-hand 
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Figure 5. A demonstration of the degeneracy between To and 7 in determining the small-scale structure of the Lya forest. The panels 
are the same as in Figure [2] except that run D15 has been replaced with run D10, which has a smaller value of 7. At z = 2.553 the Lya 
forest is sensitive to moderate ovcrdensities, where runs B15 (blue lines) and D10 (orange lines) have similar temperature. Consequently, 
the Lya forests appear very similar and have similar curvatures. At z = 4.586, the forest is more sensitive to gas near the mean density. 
A change in 7 therefore does less to alter the curvature. 



panels of Figure [6l in contrast, we plot the temperature at 
an optimal overdensity that has been empirically chosen to 
minimize the variations in the relationship between T(A) 
and log(\hi\). In this case, T(A) is a nearly one-to-one func- 
tion of log{\K,\) regardless of 7. These values can be inter- 
preted as characteristic overdensities probed the Lya forest, 
though in practice the curvature measurement will be an av- 
erage over a range of densities. The important point is that 
we can we can identify a value of A where T(A) is a simple 
function of the curvature, which allows us to reliably mea- 
sure the temperature at a density probed by the Lya forest 
without knowing the full form of the temperature-density 
relation. iMcDonald et all (|200ll 'l used a similar approach in 
their line-fitting analysis, fitting T and 7 at an overdensity 
where the errors were uncorrelated. 

The optimal overdensities for the curvature analysis will 
depend somewhat on the details of the method and on the 
quality of the data. We have selected our overdensities af- 
ter computing from the simulations using the method 
outlined in Section [4] The results are shown in Figure [7| 
We computed the optimal overdensities once using noise- 
free spectra and again after noise was added to the simula- 
tions at levels consistent with the real data. The results are 
very similar, although adding noise tends to produce a slight 
bias towards higher overdensities at z > 3. In our analysis 
we use the optimal overdensities computed with noise, and 
determine the temperatures based on the curvature values 
of the simulations with 7 ~ 1.5. Deviations from the T(A)- 
(|k|) relationship for other values of 7 are typically less than 
200-400 K, and are always less than our final la uncertainty 
in T(A). Finally, we note that even at redshifts where the 
optimal overdensity is several times the mean density the 
temperature measurement will be dominated by gas that is 



photoheated, although shock heating begins to make a mi- 
nor contribution at z < 2. 



4.4 Summary of Method 

Our method for measuring temperatures from the curvature 
of the Lya forest is thus summarized as follows: 

(i) A b-spline is fit to the data. We use an initial break 
point spacing of 50 km s~\ adding more break points as 
needed to achieve a good fit. The minimum break point 
spacing is set by the instrumental FWHM. 

(ii) The spectra are divided into 10 h^ 1 Mpc sections, and 
the flux is re-normalized by the maximum value of the fa- 
spline fit in each section. 

(iii) The curvature, k, is computed from the re- 
normalized b-spline fit. 

(iv) The path length-weighted mean absolute curvature, 
is computed from all unmasked pixels where the re- 
normalized fit falls in the range 0.1 < F R1 ° < 0.9. We group 
the data in redshift bins of Az = 0.4. 

(v) For each 10 /i Mpc section of real data, the sim- 
ulated spectra from each simulation run are smoothed to 
the same instrumental resolution and rebinned to the same 
pixel size. The (\k\) values are then computed from the 
synthetic spectra using the procedure described above after 
noise has been added and b-spline fits made. To save com- 
putational time, a grid of noise values are used (r.m.s. = 
[0.005,0.01,0.02,0.04,0.08]) and the values of (|«|) are in- 
terpolated between noise levels and simulation redshifts to 
match the data. A path length-weighted (\k\) is computed 
for each simulation using the same weighting as was used 
for the real data in each redshift bin. 
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Figure 6. Left panels: To as a function of 2o(?(|re|) for our main 
simulation runs. Colors correspond to different values of To, as 
in Figure [T] Symbols denote different values of 7 - solid circles: 
7 ~ 1.5, open circle: 7 ~ 1.3, boxes: 7 ~ 1.0, star: 7 ~ 0.7. The 
correspondence between the curvature and To becomes increas- 
ingly confused at low redshifts, where not knowing 7 can lead 
to a very wrong estimate of To. Right panels: The same, except 
that the temperatures at each redshift are measured at an opti- 
mal overdensity that has been chosen such that T(A) is a nearly 
one-to-one function of Jog(|re|), regardless of 7. 



(vi) T(A) is determined by interpolating the T(A) verses 
log{\n\) relationship in the simulations to the the value of 
log{\n\) measured from the data. The overdensity, A, is the 
optimal overdensity plotted in Figure [7] 



4.5 Impact of Thermal History 

The small-scale structure of the Lycv forest, as measured 
by the curvature or any other method, will depend on the 
thermal broadening of the absorption features as well as 
on the Hubble broadening and any turbulent broadening. 
Heating the gas will increase not only the thermal broad- 
ening but also the characteristic physical size of absorbers, 
an effect known as Jeans smoothing. The result is an addi- 
tional smoothing of the Lycv forest due to greater Hubble 
broadening across individual absorbers. In principle, there- 
fore, the small-scale structure of the Lya forest will depend 
both on the instantaneous temp erature of the gas an d on its 
integrated thermal history (e.g.. IPawlik et all 20091 '). 



K 4 



With Noise 
Without Noise 



Figure 7. The optimal overdensity at which T(A) is a nearly 
one-to-one function of the curvature. The solid line shows A (2) 
for noise-free spectra, while the dotted line shows A(jz) for spectra 
with noise properties similar to the real data. The noise ampli- 
tude increases with redshift, where it creates a mild bias in the 
curvature measurement towards higher overdensitics. 




Figure 8. Thermal histories of a subset of simulations used in 
this work. The solid lines are for runs A15-E15, with colors cor- 
responding to Figure [T] The gas temperature in these runs is set 
by photoheating in ionization equilibrium, but does not include 
an initial temperature boost due to hydrogen rcionization. The 
dotted lines shows the temperature at the mean density for runs 
T15slow and T15fast, which include heating meant to mimic an 
instantaneous hydrogen rcionization at z = 9, followed by either 
an extended (T15slow) or rapid (T15fast) He II rcionization com- 
pleting near z ~ 3. We apply our curvature method to artificial 
data drawn from these two runs in order to test the impact of the 
thermal history on measurements of the instantaneous tempera- 
ture (see Section IT5t . 



In order to test the sensitivity of our method to the 
thermal history of the IGM, we applied our analysis to a 
set of spectra drawn from simulations whose thermal his- 
tories included simple scenarios for hydrogen and helium 
reionization. The test runs, T15slow and T15fast, include 
heating meant to mimic H I reionization at z = 9, and ei- 
ther an extended He 11 reionization spanning 3 < z < 5 
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Figure 9. Temperatures recovered from the test runs using our method. Left panels are for run T15slow, right panels for T15fast. Top 
panels: Temperature at the optimal overdensity. The overdensities at which T(A) is measured are printed along the bottom of each plot. 
Bottom panesl: Temperature at the mean density inferred for the fiducial 7 ~ 1.5, which was used in both test runs. In all panels, points 
with error bars (2a) give the values we measure using artificial spectra with resolution, noise properties, and sample size similar to the 
real data. The solid lines give the actual values measured directly from the simulation. The dashed and dotted lines in the bottom panels 
give the values of To that would be recovered from noise-free HIRES and MIKE data, respectively. Differences in the thermal history of 
the test runs and the comparison simulations (cf. Figure [5} cause the increase in To to be underestimated by ~4 000K. 



(T15slow) or a rapid He 11 reionization over 3 < z < 3.5 
(TT5fast). Current estimates of the ionizing emissivity from 
QSOs suggest that He n reionization should be a patchy 
and extended process, producing a gr adual evolution in 
the volume-averaged temperature (e.g., iBolton et al.l [20091 : 
iMcQuinn et aUboOflh . Nevertheless, we include the T15fast 
model in order to demonstrate that our method is sensitive 
to rapid changes in temperature, despite the fact that the 
temperatures of the comparison models evolve slowly. For 
the sake of simplicity, the slope of the temperature-density 
relation is held to the fiducial value (7 ~ 1.5) in these runs. 
The temperature histories of the test runs are compared to 
runs A15-E15 in Figure[8] Artificial data were generated us- 
ing similar redshifts, noise levels, instrumental resolutions, 
and sample sizes as the real data. The curvature was then 
measured using the same procedure described in Section |H 
with the exception that no metal masking was performed. 

The temperatures recovered from the test simulations 



are plotted in Figure [9] The top panels shows the tempera- 
ture at the optimal overdensity, T(A), with the value of A 
at each redshift printed along the bottom of each plot. The 
bottom panels shows To, which we have computed using the 
values of 7 measured in the simulations. In all panels the re- 
covered temperatures are plotted using filled symbols, while 
the actual temperature measured directly from the simula- 
tions is shown with a solid line. We also show the values of To 
that would be recovered from noise-free spectra. The results 
at HIRES (6.7 km s" 1 ) and MIKE (13.6 km s" 1 ) resolu- 
tion are plotted using dashed and dotted lines, respectively, 
where the smoothing has been applied to both the test sim- 
ulations and the comparison runs. In general, the curvature 
method recovers the true values very well. Notably, a grad- 
ual increase in To from 2 ~ 5 to 3 is correctly recovered 
for T15slow, while a sudden jump in To at z ~ 3.5 is re- 
covered for T15fast. In both cases, however, the measured 
temperatures show some deviation from the true values. At 
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z ~ 5, T(A) is overestimated by ~2 000K, while at z ~ 3, 
T(A) is underestimated by ~3 000K for both T15slow and 
T15fast. The corresponding net rise in To is underestimated 
by ~4 000K. 

The errors in the recovered temperatures can be under- 
stood in terms of the differences in thermal history between 
the test runs and runs A15-G15, which were used to cali- 
brate the conversion from to T(A). At z ~ 5, the T15 
runs (which are identical at z > 5) and run B15 have similar 
values of To- Since the T15 runs included a heat injection 
at z ~ 9, however, the Jean smoothing in these runs will be 
substantially higher than in run B15. This should cause the 
T15 runs to have a comparatively lower overall curvature 
at z ~ 5, leading to a higher estimate of the instantaneous 
temperature. We note, however, that at z ~ 5 the impact of 
the thermal history on the measured temperature depends 
partially on the properties of the data. The recovered value 
of To in the highest-redshift bin is very near the true value 
(solid line) for noise-free HIRES data (dashed line), but in- 
creases at MIKE resolution (dotted line) or when noise is 
added (filled point, for which 90 per cent of the spectra are at 
HIRES resolution). These effects relate to properties of the 
Lya forest that are sensitive to the amount of Jeans smooth- 
ing. For a given overall mean flux, a lower amount of Jeans 
smoothing at z ~ 5 will tend to produce sharper and more 
distinct transmission peaks, which will dominate the curva- 
ture signal. With increased Jeans smoothing, more of the 
high-curvature regions for the same To will instead be small 
undulations on top of broader transmission regions. Such 
undulations are more easily obscured than distinct peaks by 
smoothing or adding noise. The curvature at z ~ 5 in the 
T15 runs is therefore lower than in run B15 when measured 
from present-quality spectra, and so the temperature in the 
test runs is over-estimated. 

The impact of thermal history is more straightforward 
at z ~ 3. Here, changes in the temperature measurement due 
to noise and lowered resolution are minimal. The T15 runs 
and run D15 have similar values of To, but because the gas in 
the T15 runs was colder than in run D15 at higher redshifts 
it has experienced less Jeans smoothing. The curvature in 
the T15 runs at z ~ 3 is therefore higher by comparison, 
which leads to an underestimate of the temperature. 

Although the instantaneous temperature is exactly re- 
covered at certain redshifts, this does necessarily mean that 
the gas has returned to hydrostatic equilibrium at these 
points. The Lya forest absorbers at a given density will tend 
to have physical sizes on order of the local Jeans length, 
Lj ~ c s tdyn, wher e e s is the sou nd speed and fd yn is the dy- 
namical timescale (|Schavell200ih . Following a heating event, 
therefore, local hydrostatic equilibrium will only be restored 
on the dynamical timescale, which at the mean density is on 
the order of a Hubble time: ta yn = l/y/Gp ~ tn A -1 . Even 
at z = 2 — 3, where the forest is probing mildly overdense re- 
gions, the thermal response time for typical absorbers will be 
significantly longer than the timescales over which hydrogen 
and helium reionization occur. In that case, the small-scale 
structure of the Lya forest is always sensitive to the thermal 
history of the gas. The redshifts where we exactly recover 
the temperatures in the test simulations are simply those 
where the integrated heat input in the T15 runs has pro- 
duced a similar amount of Jeans smoothing as in the com- 
parison run with the same instantaneous temperature. We 
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Figure 10. Raw curvature measurements. Boxes with error bars 
are for the real data. Vertical error bars are 2a, while horizontal 
error bars show the redshift range spanned by each bin. Lines 
show the curvature from the simulations using the same mix of 
resolution and noise properties present in the data at each red- 
shift. Variations in the resolution and noise cause the curvature 
in the simulations to be non-smooth functions of redshift. Line 
styles and colors correspond to Figure [T] 



will revisit the implications of this test for our temperature 
measurements from the real data in the next section. 



5 RESULTS 
5.1 Curvature 

We first present our results for the raw curvature measure- 
ments. In Figure llOl we plot log{\n\) as a function of redshift. 
The data points are plotted with 2a vertical error bars in 
redshift bins of Az = 0.4. The uncertainties were estimated 
from the simulations using sets of artificial data with sim- 
ilar redshifts, noise levels, spectral resolution, and sample 
size as the real data. In redshift bins containing large num- 
bers (A/ - > 100) of data points, these agreed very well with 
bootstrap errors generated directly from the data. Values of 
log(\n\) from the simulations are also plotted. These were 
also calculated using a similar mix of redshifts, noise lev- 
els, and spectral resolutions as the real data but from the 
full set of simulated sight lines. Variations in the properties 
of the data introduce deviations from the otherwise smooth 
evolution of the curvature with redshift (cf. Figure [T]). The 
uncertainties in our measured values of log{\K\) are typically 
< 20% of the separation between adjacent simulations with 
7 ~ 1.5, though the error is somewhat larger at z = 4.8, 
where the sample size is significantly smaller. 

As a consistency check, we have also computed log{\n\) 
as a function of the re-normalized flux level. The results for 
each redshift bin are shown in Figure [TT1 Uncertainties in the 
data points were again estimated from the simulations. The 
errors in neighboring flux bins will be highly correlated, but 
we note that the shape of £o<?(|«|) verses F mo in the sim- 
ulations is in good agreement with that of the real data. 
While we use the overall mean absolute curvature to mea- 
sure the temperatures, this agreement provides assurance 
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Figure 11. Mean absolute curvature values as a function of re-normalized flux. Boxes with error bars are for the real data. The error 
bars are 2a, though adjacent flux bins are highly correlated. Lines give values from the simulations using the same mix of resolution and 
noise properties present in the data at each redshift. In general, the observed dependence of £og{|re|) on flux is well reproduced by the 
simulations. Line styles and colors correspond to Figure [l] 



that the simulated spectra exhibit realistic curvature prop- 
erties. We note that the shape of the curvature as a function 
of flux is moderately dependent on 7. This is especially true 
at lower redshifts, where smaller values of 7 produce higher 
peaks near F R1 ° = 0.1 and a more pronounced inflection 
near F R1 ° = 0.4. Setting constraints on 7 is presently some- 
what difficult due to the correlations in the measurement 
errors and the fact that we have only run a small number 
of simulations with 7 < 1.5. We therefore leave this for a 
future work. 

5.2 Temperature at the optimal overdensity 

The primary result from this study is the evolution of the 
IGM temperature at the optimal overdensity probed by the 
Lya forest, which we infer directly from the curvature mea- 
surements. Our values for T(A) are given in Table [3] and 
plotted in Figure [T2l For each redshift bin we print the op- 
timal overdensity at which T(A) is measured along the bot- 
tom of the plot. The 2a errors are again estimated using the 
simulations, and include the errors related to uncertainties 
in Teffl The median error is 8 per cent, although this does 
not include uncertainties related to the thermal history. 

A strong trend of increasing T(A) with decreasing red- 
shift is evident. Unless the temperature-density relation is 
flat or inverted over this entire interval, which appears un- 
likely, much of this rise in temperature is probably due to the 

2 These tend to be small since we are restricting our analysis to 
the center of the flux distribution. The mean absolute curvature 
does vary with flux, as shown in Figurc im However, the shape of 
the flux PDF in the range 0.1 < F < 0.9 does not change quickly 
as Teff changes, and the flux- weighted (|k|) therefore remains rel- 
atively constant. The 2a error in T(A) related to uncertainty in 
T e ff is less than 400 K at 2 > 2.8, increasing to ~1300 K at z = 2.0. 



fact that the overdensity probed by the Lya forest data with 
our measurements increases with redshift, from A = 1.2 at 
z = 4.8 to A = 5.7 at z = 2.0. The lower density gas is adia- 
batically cooled as the voids expand, whereas the moderate 
density gas is heated as the filaments collapse. For the rel- 
atively steep temperature-density relation expected follow- 
ing hydrogen reionization, the higher density gas probed at 
lower redshifts is therefore expected to be hotter. Interpret- 
ing the T(A) results in the context of He n reionization is 
thus not straightforward, and requires a comparison to sim- 
ulations of He 11 reionization in which the full temperature- 
density relation is modeled. 

5.3 Temperature at the mean density 

Translating our T(A) measurements into temperatures 
at the mean density requires knowing the slope of the 
temperature-density relation. While we have not attempted 
to fit the full relation, we can infer what To should be for 
reasonable values of 7. We plot To in Figure [T3l under two 
different assumptions. In the first case, we use the values of 
7 measured in our fiducial simulations (A15-G15), where 7 
increases from 7 « 1.4 at z = 5 to 7 ~ 1.6 at 2 = 2. The 
exact values for 7 are given by the solid lines in the middle 
panel of Figure [T] and in Table [3] In this case, shown by 
the filled circles, To is consistent with ~8 000K at z > 4.4, 
increasing to ~12 000K at z — 2.8, and then decreasing to 
~11000K at = 2. This scenario represents the minimum To 
case, in which there is no additional heating in the voids 
that would tend to flatten the temperature-density relation. 
Even in this case, an increase in To between z — 4.4 and 
z = 2.8 is detected at the 7a level. 

In the second case, we assume a constant 7 = 1.3. 
A value of 7 in this range is sugges ted by the numerical 
simulations of iMcQuinn et all |2009l ). It corresponds to a 
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Table 3. Temperature results. Columns give the redshift interval, the path length-weighted mean redshift, the optimal overdensity, the 
total number of 10 ft. Mpc sections included, the temperature at the optimal overdensity, the fiducial value of 7, the temperature at 
the mean density for the fiducial 7, and the temperature at the mean density for 7 = 1.3. All errors are la. 
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Figure 12. IGM temperature at the optimal overdensity as a 
function of redshift. These are the primary measurements ob- 
tained by this study. The overdensity for each redshift bin is 
printed along the bottom of the plot. Vertical error bars are 2<r, 
and include statistical errors estimated from the simulations along 
with errors associated with the uncertainty in the mean flux. Hor- 
izontal error bars indicate the redshift range spanned by each bin. 



mild flattening of the temperature-density relation as ex- 
pected during an extended He 11 reionization process. To is 
largely insensitive to 7 at z > 4, where the optimal over- 
density at which we measure T(A) is already close to the 
mean density. At lower redshifts, however, To is substan- 
tially higher for 7 = 1.3, up to ~17000K at z = 2.0. Fol- 
lowing the end of He 11 reionization, 7 is expected to re- 
turn to its asymptotic value (7 ~ 1-5), and so our measure- 
ments are consistent with a decline in To at z < 2.8 even if 
the temperature-density relation is flattened at higher red- 
shifts. We emphasize that the 7 relevant to our study is 
the globally-averaged value, even though there may be sig- 
nificant scatter in the temperature-density relation. Within 
growing He ill bubbles, dense regions should ionize first 
and have time to cool b efore the voids become fully ionized 
(IFurlanetto & ohlboosD . Radiative transfer effects may also 
create additional complex ity in the temperature-density re- 
lation jBolton et alJlioOiih . Getting the temperature-density 
relation to flatten across the entire IGM simultaneously, 
however, would require a brief and spatially-coordinated 
He n reionization, which appears not to be favored on theo- 



Figure 13. Temperature at the mean density inferred from our 
T(A) measurements under different assumptions for 7. Filled cir- 
cles are for the fiducial values of 7 listed in Table [3] which are 
typically near 1.5. This is the maximum slope expected in pho- 
toionization equilibrium, and hence the filled circles represent 
minimum values for Tq. Open circles give To for 7 = 1.3. To 
for an isothermal IGM (7 = 1.0) would be given by the T(A) 
values in Figure 1121 Note that To at z > 4 is largely insensitive 
to 7, since the optimal overdensity at which we measure T(A) is 
very close to the mean density (cf. Figure [12} . The dotted line is 
a thermal asymptote of the form To oc (1 + 2) 53 that has been 
scaled to approximately match our measurements at z > 4.4. It 
shows the ex pected temperature e volution in the absence of He II 
reionization llHui fc H aiman 200j|)- 



retical grounds (e.g.. iMcQuinn et al.ll2009l ). Nevertheless, if 
the IGM is isothermal then To would be given by the values 
of T(A) in Figure [H 

Our test of the impact of thermal history in Section 14.51 
suggests that we may have somewhat underestimated the 
increase in To from z > 4 to z ~ 3. We emphasize, however, 
that since we are comparing the data to simulations with 
flat thermal histories, the position of the peak in To(z) for 
the 7 ~ 1.5 case in Figure [12] is likely to be correct. Future 
measurements of the spatial coherence of the absorbers us- 
ing pairs of QSO sight lines should allow us to disentangle 
the Jeans smoothin g from the instantaneous temperature 
(|Peeples et al.ll2010h . For now, the detection of an increase 
in To from z > 4 to z ~ 3 is robust, since the colder simula- 
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Figure 14. Comparison of our results with selected resul ts from the literat ure. The left-hand panel shows our To points for 7 ~ 1.5 and 
7 = 1.3 along with To from the wavelet analysis done by I Lidz et al. All vertical error bars in this panel are 2a. The right-hand 

panel shows our results along with those from the Voigt-profile analysis of lSchave et al. Vertical error bars for the Schaye et al. 

results are lcr. See text for discussion. 



tion runs are likely to underestimate the true level of Jeans 
smoothing, while in the the hotter runs it is likely to be over- 
estimated. Our analysis should therefore attribute too little 
of the change in curvature to changes in the instantaneous 
temperature, which would lead to an underestimate of the 
change in To with redshift. 

5.4 Comparison with literature values 

We compare our results with selected To values from the lit- 
erature in Figure [TH In the left-hand panel we overplot our 
values with the r esults of the wavelet analysis performed by 
iLidz et al.l(|2009h . The overall agreement is fair, with the ex- 
ception of z — 3.4, where their 2a lower bound significantly 
exceeds our inferred value of To for 7 = 1.3. To match their 
lower limit of 20 000 K, we would require a mildly inverted 
temperature-density relation (7 ~ 0.9). While tentative evi- 
dence evidence for an inverted temperature-density relation 
has been found in the flux probability d istribution function 
l|Becker et al.1 120071 ; iBolton et all 120081). it remains a the- 
oretical challenge to simultaneously achie ve this over the 
entir e IGM, even during He 11 reionization (|McQuinn et alj 
l2009h . We are marginally consistent with the Lidz et al. 
lower limit at z — 4.2. They note that their errors in this bin 
are strongly affected by uncertainties in the mean flux, and 
it is possible that a lower r e ff, such as the one we have used, 
would bring their To measurement into better agreement 
with ours. In the right-hand panel of Figure [T?] we overplot 
our values w i th the results of the Voigt profile analysis by 
ISchave et all (|200Ch . The agreement is again fair, although 



the Schaye et al. errors are only lcr. A careful measurement 
of 7 as a function of redshift will be required before our re- 
sults can be more directly compared to these works, which 
have both attempted to simultaneously constrain To and 7. 



6 IMPLICATIONS FOR HE II REIONIZATION 

We now discuss whether our temperature measurements re- 
quire that He 11 reionization occurs within 2 < z < 5. Follow- 
ing reionization, gas at the mean density is expected to cool 
monotonically. The gas will cool rapidly at first, then asymp- 
totically approach a thermal state where the temperature is 
set by the shape of the UVB. In photoionization equilibrium, 
the thermal asymptot e will be of the form To oc (1 + z) 0,53 
(|Hui fc Haimanll2003l ). If He 11 is reionized at z > 5, then 
this is the minimum rate at which we would expect the gas 
to be cooling over 2 < z < 5. We show an example ther- 
mal asymptote in Figure [13] which has been normalized to 
intersect our To measurements at z > 4.4. This trend falls 
well below our measured values of To at z < 4, even in the 
case of 7 ~ 1.5, which sets a lower limit for the temperature. 
The observed rise in To from z — 4.4 to z ~ 3 therefore re- 
quires an additional source of heating beyond what can be 
achieved in photoionization equilibrium. 

The expected temperature increase f or a parcel of gas 
due t o He 11 reionization is approximately (|Furlanetto fc Ohl 
|200S|) 

AT ~ 0.035 \ (E) (2) 
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where (E) is the mean energy of the ionizing photons. In the 
optically thin limit, the mean energy will be weighted by the 
ionization cross section, at oc E~ A . For an ionizing spectrum 
J v oc v~ a , this will produce (-E)thin = 54.4 (a + 2)" 1 eV , 
which for a quasar-like spectrum (a = 1.5; lTelfer et alj2002h 
gives ATthin ~ 4 000K. In the optically thick limit, all ion- 
izing photons will be absorbed and (-E)thick = 54.4 (a — 
l) _1 eV, or ATt h i ck « 30 000 K. The increase in T in Fig- 
ure [T3] can easily be accommodated within this range. We 
emphasize, however, that the temperatures we are measur- 
ing are averaged over large volumes. During an extended 
reionization process, the first regions to ionize will have time 
to cool before the overall process completes. The maximum 
increase in the globally-averaged To will therefore be smaller 
than the increased experienced by an individual parcel of 
gas. 

We leave detailed modeling of the rise in To for future 
work. Assuming that the heating of the low-density IGM at 
these redshifts is indeed dominated by He II photoheating, 
however, we can already formulate two basic conclusions: 
(i) that He II reionization occurs at z < 5, and (ii) that the 
process is extended, beginning at z > 4.4 and potentially 
extending to z ~ 3, consistent with the end of He II reion- 
ization inferred from the evolution in the He II mean opacity 
(Dix on fc Furlanettol |2009| ; iMcQuinnl 120091 . and references 
therein). The redshift at which To(z) begins to decline, as 
well as the peak value of To will depend primarily on the 
evolution of 7 with redshift. 

Finally, we note that our values of T> at 2 ~ 4 — 5 are 
consistent with a broad range of hydrogen reionization sce- 
narios. At z > 6, To is expected to decrease rapidly following 
H 1 reionization due to the combi ned effects of adiabatic 
expa n sion and Compton cooling (|Miralda-Escude fc Reesl 
1 1994 iHui fc Gnedmfll997l '). Even for H I reionization at 
z > 10, however, the asympt otic value of Tp at z ~ 4 — 5 
is expected to be ~8 000K (|Hui fc Haimanl |2003h . which 
is consistent with our measurements. This temperature as- 
sumes a quasar-like background with a — 1.5, while a stel- 
lar spectrum with a = 3.0 would lower To by ~20 per cent 
(|Hui fc Haimanl l2003l ). On the other hand, our estimates 
of the temperature at z > 4 may be too high by roughly 
the same amount for the reasons described in Section 14.51 
Given the uncertainties in both the measurements and the 
thermal models, therefore, we conclude that our temper- 
ature measurements at z > 4 do not obviously require a 
recent (z < 10) temperature boost f rom H I reionization . 
This in in contrast to earlier works (|Theuns et al 1 l2002al : 
IHui fc Haimanll2003h based on temperatures estimates that 
were considerably higher than the ones measured here. We 
note that the rapid cooling expected after reionization com- 
pletes means that our measurements are also consistent with 
H 1 reionization at z < 10 for a wide vari ety of reioniza- 
tion redshifts and initial temp e rature boosts (IHui fc Haimanl 
120031 ; iFurlanetto fc Ohl 12008 ). iBolton et al.l ([ioifF recently 
demonstrated that local constraints on the timing of hy- 
drogen reionization can also be placed from temperature 
measurements in the near-zones of z ~ 6 QSOs. Combin- 
ing near-zone measurements with our results from the forest 
should in the future provide a powerful means of studying 
the history hydrogen reionization. 



7 SUMMARY 

We have produced new measurements of the IGM temper- 
ature over 2.0 < z < 4.8. These are the most sensitive such 
measurements to date, and the first in the general IGM at 
z > 4.5. Our analysis uses a new statistic, the curvature, to 
characterize the small-scale structure of the Lya forest in 
a manner that is highly sensitive to the IGM temperature 
while allowing systematic errors to be effectively minimized. 
We further focus on measuring the temperature at an opti- 
mal overdensity for each redshift, where the temperature is 
nearly a one-to-one function of the mean absolute curvature. 
The temperature at the mean density is then inferred from 
our T(A) measurements for different values of the slope for 
the temperature-density relation. Our main conclusions are 
as follows: 

• The temperature at the optimal overdensity increases 
from T(A) « 9 000K at z = 4.4 - 4.8, where A = 1.2 - 1.4, 
to T(A) « 30 000 K at z = 2.0, where A = 5.7. This trend is 
most likely driven in large part by the fact that A increases 
with decreasing redshift, since the temperature of the IGM is 
generally expected to increase towards higher densities due 
to adiabatic heating in the collapsing filaments and cooling 
in the expanding voids. The T(A) measurements are the 
primary result of this study, and can be directly compared 
with theoretical models of He 11 reionization. 

• We use the maximum slope for the temperature-density 
relation expected in the absence of He 11 reionization (7 ~ 
1.4 — 1.6, depending on redshift) to set a minimum value 
of To as a function of redshift. We also consider the case of 
7 = 1.3, motivated by re cent He 11 reionization simulations 
bv lMcQ uinn et al. I(|2009h. In either case (and for any smaller 
values of 7), we detect an increase in To from z = 4.4 to 
z = 2.8. The minimum increase is ATo « 4 000K, although 
the increase may be larger if 7 decreases over this interval, 
or if the differing amounts of Jeans smoothing in our sim- 
ulations have caused us to underestimate the the change in 
temperature. 

• The temperature evolution from z ~ 4.4 to 2.8 is 
not consistent with the asymptotic temperature evolution 
of a photo-ionized IGM expected long after reionization is 
completed, where Tp decreases monotonically with redshift 
(|Hui fc Gnedinlll997T l. The observed increase in Tp over this 
interval requires a substantial injection of additional energy, 
which we attribute to photoionization heating during an ex- 
tended period of He 11 reionization. 

• A peak in To(z) at z ~ 2.8 is inferred for 7 ~ 1.5. This 
would be consistent with an end to He 11 reionization at 
z ~ 2.7, as inferred from He 11 Lya opacity measurements. 
However, the location of the peak can change if 7 evolves 
with redshift. Our measurements are consistent with a falloff 
in Tp at z < 2.8, particularly since 7 is expected to return 
to its asymptotic value (7 ~ 1.5) following the end of He 11 
reionization. 

• The value of Tp we infer from our T(A) measurements 
over z — 4.4 — 4.8 is largely independent of 7, since we 
are already measuring the temperature close to the mean 
density. At these redshifts, To « 8 000K, which is consistent 
with H 1 reionization having completed well before 2 ~ 6. 

Two key challenges remain to fully characterizing the 
thermal evolution of the IGM over this redshift range, at 
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least in a volume-averaged sensed. The most pressing issue 
is to establish the shape of the temperature-density relation 
and its evolution with redshift. This will clarify the evolu- 
tion of To(z) at z < 4, and allow us to identify the end of 
He II reionization based on the redshift at which To peaks. 
The second challenge is to determine the amount of Jeans 
smoothing in the IGM, which can possibly b e constrained 
using lines of sight towards pairs of QSOs (jPeeples et alj 
bold ). We leave both of these tasks for the future, but note 
that the measurements presented here provide a basis for 
making further progress. For example, our T(A) values may 
be combined with Voigt profile measurements to determine 
7 over z ~ 2 — 4. 

Ultimately one would like to characterize not only the 
evolution of the volume-averaged thermal state of the IGM 
but also the evolution in its scatter. A key characteristic 
of He II reionization should be that growing He ill bub- 
bles around QSOs introdu ce spatial fluctuations in both ion- 
ization and temperature (ISokasian et a l. 2002; Glcsc r et al 



20051; lLai et al.l 1200(1 IPaschos et alj|2007l ; iFurlanetto fc Oh 



20081 ). eve n if these are partially damped by radiative trans- 



fer effects (|McQuinn et al.||2009h . Establishing the volume- 
averaged evolution of IGM temperatures provides a baseline 
against which to measure these fluctuations. The timing and 
amplitude of the increase in To, the changes in 7 and the 
evolution of the scatter will each be essential ingredients for 
understanding the role of He 11 photoheating in the ther- 
mal balance of the IGM, and for formulating a more robust 
picture of when and how He 11 reionization occurred. 
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2 x 256 3 . In that case, further increasing the mass resolution 
by a factor of eight increased Ioq(\k\) by only 0.014, which at 
z — 4.915 would correspond to a change in temperature of 
roughly 500 K. While we do not apply a mass resolution cor- 
rection to our results, therefore, we can be reasonably con- 
fident that the gains from going to higher resolution would 
be small, and that they would only be significant our high- 
est redshift bins where the correction is likely to be smaller 
than the errors. 



APPENDIX A: NUMERICAL CONVERGENCE 
TESTS 

Capturing the small-scale structure in the high-redshift Lya 
forest poses a significant challenge for co s molog ical simu- 
lations. As discussed by iBolton fc Becker] (|2009h . most of 
the transmission in the forest at z > 5 comes from voids, 
which are less well resolved in SPH simulations than the 
filaments that dominate the characteristics of the forest at 
z ~ 2 — 3. Ideally, one would use both a large simulation box 
(I > 40 h' 1 Mpc), and high resolution (M gas < 10 5 h' 1 M ). 
Computational constraints, however, required us to focus on 
obtaining the necessary mass resolution, which was found to 
be a more stringent requirement for the curvature. To check 
our convergence with mass resolution and box size, we ran 
a set of comparison simulations (runs Rl - R4). In these 
runs we used the same thermal history as in run C15, but 
varied the box size and mass resolution. Compared to run 
C15, which uses L = 10 h~ x Mpc, n — 2x 512 3 particles and 
M gas = 9.2 x 10 4 h' 1 M , the mass resolution was lowered 
by up to a factor of 64, and the box size was increased by up 
to a factor of four in length. The optical depths in each sim- 
ulation were rescaled such that the mean fluxes were equal 
after the spectra were re-normalized in 10/i _1 Mpc sections 
(Section l4.2.3p . The curvature was then computed using the 
procedure outlined in Section[4]4] with the modification that 
noise-free spectra were used and no b-spline fits were per- 
formed. 

The convergence results are shown in Figure IAT1 where 
we plot the change in log{\n\) with respect to run C15. We 
are nearly converged with box size at all redshifts. At a fixed 
resolution, log(\n\) always decreases by less than 0.03 when 
increasing the box size from 10 to 40 h Mpc. For the curva- 
ture values measured in the data, this corresponds to a tem- 
perature difference of less than 900 K. A box size correction 
calculated from these results was applied to the simulation 
curvature values when fitting the temperatures. At z < 4 we 
are well converged with mass resolution, although the con- 
vergence is less clear at higher redshifts. At z = 4.915, the 
increase in log{\n\) when increasing the number of particles 
from 2 x 256 3 to 2 x 512 3 is 0.058. We can estimate the gains 
from an additional factor of eight in mass resolution by com- 
paring to z — 3.991, which has a similar change in log(\n\) 
when increasing the number of particles from 2 x 128 3 to 
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Figure Al. Convergence of the mean absolute curvature with box size and resolution. Curvature values were computed from spectra 
without noise, using the same method as was applied to the real data. The change in (\k\) is measured with respect to the highest- 
resolution case. We are only marginally converged with mass resolution at the high redshift end, but the likely change in (\k\) gained 
from further increasing the number of particles would only correspond to a small change in the measured temperature (see text). We do 
apply a box size correction to the curvature, although it is small. 
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